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ABSTRACT 

We  have  a  system  of  N  sendox  machines  (service  centers) 
which  can  transmit  or  receive  (but  not  simultaneously)  photo- 
copies to  and  from  each  other.   The  stochastic  process 

{X„(t),t  >  0},  where  X„(t)  is  the  number  of  messages  under- 
M       —  M 

going  transmission  in  the  system  at  epoch  t,  is  selected  to 
be  the  measure  of  system  congestion.   Under  the  assumptions 
of  exponential  inter-arrival  and  transmission  times  and  that 
messages  not  transmitted  immediately  are  lost,  X„(t)  is  a 
birth  and  death  process.   We  call  the  model  under  these  as- 
sumptions the  "loss"  model.   It  is  shown  that  X„(t)  can  be 
approximated  by  a  diffusion  process.   The  steady  state  re- 
sults of  the  diffusion  model  can  be  calculated  with  ease  and 
are  in  close  agreement  with  those  of  the  birth  and  death 
model  and  the  simulation  "loss"  model.   A  more  realistic  "no 
loss"  model  is  also  developed  where  messages  that  are  not 
transmitted  immediately  are  allowed  to  queue  up.   It  is  shown, 
however,  that  this  "no  loss"  model  can  be  quite  accurately 
approximated  by  an  appropriate  "loss"  model. 


TABLE  OF  CONTENTS 

I.  INTRODUCTION 10 

II.  DESCRIPTION  OF  THE  "LOSS"  MODEL 12 

III.  BIRTH  AND  DEATH  "LOSS"  MODEL  14 

IV.  DIFFUSION  "LOSS"  MODEL  16 

V.  SIMULATION  "LOSS"  MODEL 25 

A.  SIMULATION  PROGRAM  I 25 

B.  SIMULATION  PROGRAM  II 26 

C.  DISCUSSION  OF  OUTPUT  OF  PROGRAMS  I  Si  II 27 

VI.  SIMULATION  "NO  LOSS"  MODEL  29 

VII.  COMPARISON  OF  RESULTS 31 

VIII.  APPROXIMATION  OF  "NO  LOSS"  MODEL  BY  "LOSS"  MODEL-  39 

IX.  CONCLUSION 44 

APPENDIX  A  Derivation  of  Characteristic  Function  of 

S(t) 45 

APPENDIX  B  Derivation  of  Stochastic  Differential 

Equation  of  S(t)  51 

APPENDIX  C   Derivation  of  Steady  State  Mean  and 

Variance  of  X.,(t)  54 

M 

APPENDIX  D   Flowcharts  and  Listings  of  Computer  Programs  59 

APPENDIX  E   Calculation  of  A'  82 

LIST  OF  REFERENCES 84 

INITIAL  DISTRIBUTION  LIST  85 


LIST  OF  TABLES 


Table 


I.  Steady  state  mean  and  variance  of  X,(t)  with  p 
constant 32 

II.  Steady  state  mean  and  variance  of  X M(t)  with  N 
constant 34 

III.  Steady  state  mean  of  X„(t)  with  X   and  u  constant  -  35 

IV.  Steady  state  variance  of  X„(t)  with  X    and  u 
constant 36 

V.  Steady  state  cumulative  probabilities  for  the 

three  "loss"  models 37 

VI.  Comparison  of  "no  loss"  model  with  its  "equivalent" 
"loss"  model  — 41 

VII.  Steady  state  cumulative  probabilities  for  the 
"loss"  and  "no  loss"  models 42 


6 


LIST  OF  FIGURES 


Figure 

1.  Plot  of  p  vs.  x  /M 21 

2.  Plot  of  p  vs.  a2/M 22 

3.  Plot  of  steady  state  cumulative  probabilities  vs. 
X„(t)  for  the  three  "loss"  models 38 

4.  Plot  of  steady  state  cumulative  probabilities  vs. 
X.,(t)  for  the  "loss"  and  "no  loss"  models 43 

5.  Plot  of  x  vs.  B(x) 57 


NOTATIONS 

N      total  number  of  machines  (service  centers)  in  the  sys- 
tem.  (Assumed  to  be  an  even  number.) 

M=N/2   maximum  number  of  messages  that  can  be  transmitted  by 
the  system. 

A.      arrival  rate  of  messages  at  each  machine. 

X  _„    effective  arrival  rate  of  messages  at  each  machine  in 
the  "loss  model. 

u      departure  rate  of  a  message  undergoing  transmission. 

P 


y 

2A 


X„(t)  number  of  messages  undergoing  transmission  at  epoch  t 

S.,(t)  stochastic  element  associated  with  X,.(t). 
M  M 

y(t)  deterministic  component  of  X„(t). 

<J>w(t)  characteristic  function  of  X„(t). 

ik,(t)  characteristic  function  of  SM(t). 

x  steady  state  mean  of  X.,(t). 
s  J  M 

a2  steady  state  variance  of  X„(t). 

P.'s  steady  state  probabilities  of  X„(t). 
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I.   INTRODUCTION 

The  main  purpose  of  this  thesis  is  to  construct  and  com- 
pare various  models  which  can  be  used  to  analyze  the  degree 
of  congestion  in  a  system  of  N  sendox  facsimile  machines. 
Sendox  machines  can  both  transmit  and  receive  (but  not  simul- 
taneously) photocopies  using  standard  telephone  lines.   The 
original  problem  arose  out  of  a  proposal  to  install  a  number 
of  these  machines  in  central  London  Mintech  buildings  for 
transmission  of  thin  "immediate"  messages  between  buildings. 
This  problem  was  first  studied  by  Coleman  [Ref .  1]  who  in- 
vestigated the  expected  system  congestion,  expected  cost  per 
message  and  a  possible  priority  system. 

One  of  the  models  in  his  study  considered  each  machine 
individually  as  a  single  server  system.   It  was  assumed  that 
only  one  machine  was  located  in  each  building.   Arrivals  at 
each  of  the  N  machines  occur  from  two  sources:  1)  requests 
to  transmit  messages  to  other  machines  and  2)  incoming  mes- 
sages from  the  other  N-l  machines.   These  two  arrival  streams 
were  merged  into  one  and  the  standard  queuing  formula  for  a 
single  server  queue  was  applied  to  obtain  the  expected  wait- 
ing time  of  messages  at  an  individual  machine.   However,  a 
subsequent  simulation  study  showed  wide  discrepancy  between 
the  theoretical  and  simulated  expected  waiting  times.   This 
occurred  because  the  single  server  model  did  not  take  into 
account  the  "interference'  between  machines.   The  model  as- 
sumed that  when  a  message  is  transmitted  from  one  machine, 
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its  destination  (which  could  be  any  one  of  the  other  N-l  ma- 
chines) is  always  available  to  receive  the  message.   In 
reality,  its  destination  may  be  engaged  in  transmitting  or 
receiving  another  message.   The  effect  of  such  "interference" 
tends  to  increase  the  waiting  time  as  compared  to  the  result 
obtained  from  the  single  server  queuing  model.   Any  reason- 
able theoretical  model,  therefore,  has  to  consider  all  N  ma- 
chines simultaneously.   It  is  very  difficult,  however,  to 
build  an  analytical  model  that  can  give  a  reasonable  esti- 
mate of  waiting  time  because  it  has  to  consider  N  queues  all 
at  once. 
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II.   DESCRIPTION  OF  THE  "LOSS"  MODEL 

The  problem  becomes  much  more  tractable  when  we  only  con- 
sider modelling  the  system  size  i.e.,  the  number  of  messages 
undergoing  transmission.   This  number  varies  from  0  to  M 
where  M=[N/2],  the  largest  integer  less  than  or  equal  to  N/2. 
Let  X,,(t)  be  defined  as  the  number  of  messages  undergoing 

M 

transmission  at  epoch  t.   The  stochastic  process  {X„(t), 
t  >_  0}  will  be  our  measure  of  system  congestion. 

A  number  of  assumptions  are  made  to  facilitate  construc- 
tion of  theoretical  models  for  the  process  XM(t).   We  assume 
a  Poisson  arrival  of  messages  at  each  machine.   The  arrival 
rate  of  messages  requesting  transmission  at  each  machine  is 
the  same  and  is  denoted  by  A.   Moreover,  we  assume  that  each 
message  undergoing  transmission  is  finishing  at  an  exponen- 
tial rate  of  y.   Furthermore  we  define  p  to  be  -My  and  hence 
the  departure  rate  (y)  of  a  message  is  also  equal  to  2pA. 
Finally,  we  assume  that  if  a  message  fails  to  be  transmitted 
immediately  upon  arrival  (due  to  either  the  receiving  or 
sending  machine  being  busy),  then  the  message  is  lost  to  the 
system.   For  this  reason,  we  call  the  model  based  on  the 
above  assumptions  the  "loss"  model. 

Two  theoretical  and  two  simulation  models  for  the  proc- 
ess X.,(t)  will  be  described  and  discussed  in  the  following 
sections.   We  are  mainly  interested  in  steady  state  behavior 
of  the  process  X„(t).   In  the  first  theoretical  model,  X,,(t) 

M  M 

is  shown  to  be  a  birth  and  death  process  (this  was  done  in 


12 


[Ref.  1])  whose  steady  state  probabilities  can  be  derived 
and  calculated  numerically .   Since  Iglehart  [Ref.  3]  had 
proved  the  convergence  of  certain  birth  and  death  processes 
to  diffusion  processes,  we  were  led  to  consider  a  diffusion 
approximation  for  the  process  XM(t),  when  M  is  large.   This 
was  the  approach  taken  in  the  second  theoretical  model  where 
it  was  shown  that  X„(t)  can  be  approximated  by  a  non-station- 
ary Ornstein-Uhlenbeck  process.   Finally,  a  simulation  model 
was  built  to  validate  the  theoretical  models.   The  steady 
state  results  of  the  theoretical  and  the  simulation  models 
compared  very  favorably  with  one  another. 

A  second  simulation  model  was  developed  for  the  case  in 
which  messages  not  transmitted  immediately,  queue  up  until 
they  are  finally  transmitted.   We  call  this  model  the  "no 
loss"  model.   We  are  unable  to  construct  a  theoretical  model 
for  this  case.   However  we  will  show  that  this  "no  loss" 
model  can  be  quite  accurately  approximated  by  the  diffusion 
"loss"  model. 
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III.   BIRTH  AND  DEATH  "LOSS"  MODEL 

In  an  unpublished  paper  [Ref.  1],  Dr.  Rodney  Coleman 
modelled  the  number  of  messages  undergoing  transmission, 
X.,(t),  as  a  birth  and  death  process.   For  the  sake  of  sim- 
plicity,  we  assume  that  N,  the  total  number  of  machines  in 
the  system  is  even  and  hence  M,  the  maximum  number  of  mes- 
sages that  can  be  carried  by  the  system  is  given  by  N/2. 
When  the  process  X.,(t)  is  in  state  j,  the  conditional  proba- 
bility that  X„(t)  will  visit  state  j+1  in  the  next  interval 
of  length  At  is  X  .At  +  o(At)  and  the  conditional  probability 

«J 

that  X..(t)  will  visit  state  j-1  is  u  .At  +  o(At),  where  X. 
Mv  3  3 

and  u  .  are  the  infinitesimal  birth  and  death  rates  respec- 
3 

tively  [Ref.  6].   Messages  can  only  enter  the  system  at  the 
(2M-2j)  free  machines  at  a  rate  of  (2M-2j)X.   When  a  message 

arrives  at  a  free  machine,  its  probability  of  finding  an 

.,  .,       .  .       .  .    .   2M-2J-1   ..     „ 
available  receiving  machine  is  — n.j°  = — ,  therefore. 

2M-1  ' 

l   _  (2M-2j)(2M-2j-l)X 
j        (2M-1)  U,1J 

Since  there  are  j  messages,  each  finishing  transmission  at 
the  rate  of  y,  and  thus 

U-i  =  3V    =   2jpX  (3.2) 
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Let  P.  denote  the  steady  state  probability  of  X  (t)  being  in 
state  j  .   Then  Coleman  [Ref.  1]  showed  that 

k 


„  k!(2M-2k)! k=0,l,...M 

k    M      k  '  ' 

V    JL 


k=0  k!(2M-2k)! 


where  y  =  (-^j-t)  ~   - 


2M-ly  y    2p(2M-l) 

However,  since  the  caluclation  for  the  P.  's  are  carried 

k 

out  on  a  digital  computer,  a  more  efficient  method  computa- 
tionally is  to  use  the  following  recursive  relationships 
[Ref.  6]. 

ykPk  =  Xk-lPk-l       k=l,2,...M  (3.3) 

M 
kI0Pk  -  1       "  (3.4) 

The  computer  program  which  calculates  the  steady  state  prob- 
abilities and  also  the  mean  and  variance  of  X„(t)  at  steady 
state  is  listed  in  Appendix  D. 
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IV.   DIFFUSION  "LOSS"  MODEL 


We  follow  Gaver ' s  [Ref.  4]  approach  by  approximating  the 

process  X„(t)  with  a  deterministic  component  y(t)  and  a 

stochastic  element  or  noise  S.,(t).   Thus  X„("t)  can  be  written 

M  M 

as: 

XM(t)  =  My(t)  +  /M  SM(t)  (4.1) 

We  show  that  as  we  let  M  go  to  infinity,  SM(t)  converges 
to  a  non-stationary  Ornstein-Uhlenbeck  process,  S(t)  [Ref. 
2] .   Hence  for  large  M,  we  can  use  the  following  formula  to 
approximate  X„(t): 

XM(t)  -  My(t)  +  M   S(t)  (4.2) 

Because  of  our  assumption  of  Poisson  arrivals,  a  message 
will  arrive  at  a  machine  with  probability  Adt  +  o(dt)  in  the 
time  interval  (t,t+dt).   During  this  time  interval,  the 
probability  of  two  or  more  messages  arriving  at  the  same  ma- 
chine is  o(dt).   At  epoch  t  the  number  of  messages  under- 
going transmission  is  XM(t)  and  hence  the  number  of  unoccupied 
machines  is  given  by  N-2XM(t)  or  2M-2X  (t).   In  the  time  in- 
terval (t,t+dt),  effective  arrival  of  messages  can  occur 
only  at  any  of  the  2M-2X.,(t)  unoccupied  machines.   Given  that 
a  message  has  arrived  at  an  available  machine,  its  probabil- 
ity of  finding  the  receiving  machine  not  in  use  is  now 
(2M-2XM(t)-l)/(2M-l) .   Therefore, 

Prob{XM(t+dt)=XM(t)+l|given  XM(t)} 

(2M-2X„(t))(2M-2X„(t)-l)Adt  (4.3) 
2M=1 +  °(dt) 
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Similarly,  under  the  assumption  of  exponential  transmis- 
sion times,  each  message  undergoing  transmission  is  finish- 
ing at  the  rate  of  u  or  2pX.   At  epoch  t,  messages  will  be 
finishing  at  a  rate  of  2XM(t)pX.   Therefore  during  the  time 
interval  (t,t+dt)  the  probability  of  a  message  completing 
its  transmission  is  2XM(t)pXdt  +  o(dt)  and  the  probability 
of  two  or  more  messages  completing  their  transmission  is 
o(dt).   Therefore, 

Prob  {XM(t+dt)  =  XM(t)-l|  given  X^t)} 

(4.4) 
=  2X  pXdt  +  o(dt) 

Finally, 

Prob  (XM(t+dt)  =  XM(t)|  given  X^t)} 

(2M-2XM(t))(2M-2XM(t)-l)  (4-5) 

=  1 M  2M_X Xdt  -  2XM(t)pAdt  +  o(dt) 

In  Appendix  A,  equations  (4.1)  through  (4.5)  are  used  to 
derive  the  partial  differential  equation 


U  =  -2X{2  +  p-2y(t)}e||  -  2X{y2(t)-y(t)(2-p)  +  l}-|^ 


(4.6) 

Here  4>(6,t)  is  the  characteristic  function  of  S(t),  the  lim- 
it of  the  noise  process  S„(t)  as  M  ■*  °°.   The  deterministic 
process  y(t)  must  satisfy  the  ordinary  differential  equation 

y*(t)  -  2X{y2(t)  -  (2+p)y(t)  -  1}  -  0  (4.7) 

in  order  that  the  partial  differential  equation  (4.6)  be 
valid.   The  partial  differential  equation  (4.6)  is  recognized 
as  the  transformed  version  of  the  forward  differential 
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equation  for  a  non-stationary  Ornstein-Uhlenbeck  process. 
It  is  known  [Ref.  2]  that  the  probability  density  function 
of  the  Ornstein-Uhlenbeck  process  S(t)  must  satisfy  the 
Kolmogorov  Forward  Differential  Equation 

|f  -  h   ajp-  {o(s,t)p}  -  f^  {B(s,t)p}  (4.8) 

where  3(s,t)  and  a(s,t)  denote  the  infinitesimal  mean  and 
variance  respectively,  conditioned  on  S(t)=s. 

3(s,t)  =  -2A{2+p-2y(t)}s  (4.9) 

a(s,t)  =  2A{y2(t)-(2-p)y(t)+l}  (4.10) 

The  first  step  toward  determining  the  mean  and  variance 
of  the  process  S(t)  involves  solving  the  first  order  differ- 
ential equation  (4.7)  to  obtain  y(t)  explicitly,  assuming 
the  initial  condition  that  y(o)=0.   The  solution  is  derived 
in  Appendix  A  and  is 

a-4Xwt 
-4Awt 


y(t)  =  v-w  {1  +  ke }  (4.11) 


1  -  ke 

where   v  =  1  +  ip 

w2  =  (1  +  ip)2  -  1 

,    w-v 

k  =  — — 
w+v 

Since  the  Ornstein-Uhlenbeck  process  is  Gaussian,  we  may 
differentiate  the  Gaussian  characteristic  function 
exp  [  i8a(t  )--£02b(t  )  ]  with  respect  to  6  and  t,  substitute  \p , 
7T,  ■jrj;  back  into  equation  (4.6)  and  equate  the  coefficients 
of  i0  and  82  in  order  to  obtain  two  first  order  differential 
equations  for  a(t)  and  b(t)  which  denote  the  mean  and 
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variance  of  the  process  S(t)  respectively.   Assuming  that 
S(0)=0,  we  can  show  that 


E[S(t)]  -  0        t  >  0  (4.12) 


and  since 


XM(t)  -  My(t)  +  /M  S(t) 

we  obtain 

E[XM(t)]  -  My(t)    t  >  0  (4.13) 

where  y(t)  is  given  by  equation  (4.11) 

The  solution  of  the  first  order  differential  equation  for 

the  Var  {S(t)}  involves  integrals  that  are  quite  intractable 

and  we  are  only  able  to  obtain  an  approximation  for  large  t. 

y2  -  y  (2-p)+l 
Var  {S(t»     2(2+p-2ye) *  large  (4'14) 


lim 


where   ye  =  £^  y(t)  =  1  +  £p   -    /(l+ip)2-l  (4.15) 

Since  Var  {S(t)}  -  tt  Var  {X„(t)},  we  obtain  the  following 

M       M 

expressions  for  steady  state  mean  and  variance  of  XM(t),  de- 
noted by  x   and  a2  respectively, 
s       s 


xs  =  EtXM(t)]  =  Mye      t    large 
o2s   -  Var  {XM(t)} 


(4.16) 


y2-y  (2-p)+l 
=  M  t2!2+p-2ye)   ]    t   large  (4-17) 

Detailed  derivation  of  equations  (4.12)  and  (4.14)  is  given 

in  Appendix  C.   When  t  is  large,  X,,(t)  is  normally  distri- 

M 

buted  with  mean  and  variance  given  by  equations  (4.16)  and 
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(4. 17)respectively .   Plots  of  x  /M  and  a2/M  against  P  are 

s        s 

given  in  Figures  1  and  2  respectively.  It  is  seen  in  Appen- 
dix C  that  y(t)  reaches  its  steady  state  value  of  y  consid- 
erably-faster than  Var  {S(t)},  thus  the  steady  state  mean  of 
X.,(t)  is  attained  much  faster  than  the  steady  state  variance 

M 

of  XM(t). 

The  steady  state  probabilities  for  X,.(t)  can  be  easily 

M 

obtained  from  the  standard  normal  table  after  x   and  a2  are 

s      s 

calculated  using  equations  (4.16)  and  (4.17).   However,  since 

we  are  using  a  continuous  (normal)  probability  distribution 

to  approximate  a  discrete  distribution,  we  have  to  make  the 

following  "continuity  corrections" 

(x-x    )2 

1  o 

2  1  9a2 

P      »      /     — - —   e      zos  dx  (4.18) 

O  — oo  nr — 

/2tto 

S  /  >2 

(x-x   )z 

1   *         1 2a1 

P.    -    ./,   — - —  e        ^s  dx        i=l,2,...M-l  (4.19) 

1         1_2    /2Ta 

S  /  s  2 

(x-x   V 

oo  7^ 

PM  =      /     — - —     e      2as  dx  (4.20) 

M    2    /2tto_ 

S 


Another  way  to  derive  the  parameters  for  the  Ornstein- 
Uhlenbeck  process  S(t)  is  by  defining  the  stochastic  differ- 
ential equation  for  the  process  XM(t). 

dXM(t)  =  B  (x,t)dt  +  Z(t)  AxM(x,t)dt  (4.21) 

where  Z(t)  is  a  purely  random  Gaussian  process  with  zero 

mean  and  unit  variance  and  8.,(x,t)dt  and  u.,(x,t)dt  are  the 

mean  and  variance  of  the  increment  dX..(t)  in  a  small  time 

M 
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interval  (t,t+dt).   Recall  that  we  have  derived  the  follow- 
ing probabilities  earlier 

Prob  (dX„(t)  =  l|  given  X.,(t)} 

M  M 

{2U-2X    (t)}{2M-2X  (t)-l} 

= *w=i M Adt  +  °(dt) 

Prob  (dXM(t)  =  -l|  given  XM(t)} 
=  2XM(t)pXdt  +  o(dt) 

Prob  (dXM(t)  =  0|  given  *M(t)} 

{ 2M-2X  ( t ) } { 2M-2X  ( t  )-l } 

=  1 5L_ B xdt  -  2X„(t)pXdt  +  o(dt) 

2M-1  M 

Therefore  it  is  reasonable  to  define  B,,(x,t)dt  and  a„(x,t)dt 

M  M 

as  follows: 


3M(x,t)dt  =  E[dXM(t)] 


{2M-2X  (t)}(2M-2X  (t)-l) 

L_ M Xdt-2XM(t)pAdt 


(4.22) 


aM(x,t)dt  =  Var  [dXM(t)] 


(4.23) 


{2M-2X  (t)}{2M-2X  (t)-l> 
=  2M_! " Xdt  +  2XM(t)pAdt 

Substituting  SM(t)  as  given  by  equation  (4.1)  in  place  of 
X.,(t),  we  derive  a  stochastic  differential  equation  for 
S.,(t)  in  Appendix  B.   Then  by  letting  M  ■*•  °°  and  requiring 
the  first  order  differential  equation  (4.7)  for  the  deter- 
ministic process  y(t)  to  hold,  we  obtain  the  stochastic 
differential  equation 
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dS(t)  =  2A{2y(t)-2-p}S(t)dt 


+  Z(t)  /2A{y2(t)-y(t)(2-p)+l}dt  (4.24) 

This  equation  is  the  stochastic  differential  equation  for  an 
Ornstein-Uhlenbeck  process  which  has  the  same  infinitesimal 
mean  and  variance  as  those  given  in  equations  (4.9)  and  (4.10) 
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V.   SIMULATION  "LOSS"  MODEL 

A.   SIMULATION  PROGRAM  I 

Two  methods  of  simulations  are  presented  here  for  the 
"loss"  model.   Simulation  Program  I  makes  use  of  the  obvious 
approach   of  duplicating  the  operation  of  the  actual  system. 
Messages  arrive  at  each  of  the  N  machines  at  a  Poisson  rate 
X.   These  messages  are  either  lost  or  undergo  transmissions 
which  are  being  completed  at  an  individual  rate  of  u .   A 
brief  description  of  the  program  logic  is  given  in  the  next 
paragraph. 

Simulation  Program  I  is  written  in  FORTRAN.   It  operates 
simply  by  searching  for  the  next  event  to  occur.   In  this 
model  an  event  is  either  the  arrival  of  a  message  at  a  ma- 
chine or  the  completion  of  a  message  undergoing  transmission, 
Each  event  is  assigned  a  clock  time  which  denotes  the  time 
of  occurrence  of  that  event.   The  next  event  is  that  whose 
clock  time  is  the  minimum  of  the  clock  times  of  all  the  fu- 
ture events.   Having  found  this  event,  the  program  performs 
a  series  of  tests  to  determine  whether  the  conditions  in  the 
system  are  such  that  the  potential  event  can  take  place.   If 
so,  the  program  moves  to  execute  this  event.   For  example 
the  event  may  be  the  arrival  of  a  message  at  machine  i,  then 
the  necessary  test  is  to  check  whether  the  sending  machine 
i  and  a  randomly  selected  receiving  machine  j ,  jfx    are  both 
available.   If  so,  the  number  of  messages  undergoing  trans- 
mission increases  by  one.   The  basic  loop  is  repeated  as 
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many  times  as  necessary  to  complete  the  entire  simulation. 
The  flow  chart  and  listing  of  this  program  are  given  in  Ap- 
pendix D. 

B.   SIMULATION  PROGRAM  II 

Another  way  to  develop  a  simulation  program  for  the  "loss" 

model  is  to  make  use  of  the  fact  that  X.,(t)  is  a  birth  and 

M 

death  process.   We  know  the  process  X..(t)  stays  in  state  i 
for  an  amount  of  time  that  is  exponential  with  rate  A.+u. 

and  given  that  it  makes  a  transition  from  state  i,  it  goes 

Vi 

to  either  state  i-1  or  i+1  with  probabilities  -r — 7 —  and 

■r — - —  respectively  where  A.  and  u.  are  as  defined  in  section 

III. 

Instead  of  monitoring  the  flow  of  messages  through  the 

system  of  N  machines,  Program  II  works  by  keeping  track  of 

the  length  of  time  that  the  process  X„(t)  spends  in  each 

state  and  when  XM(t)  changes  state.   This  program  is  also 

written  in  FORTRAN.   It  starts  off  with  X  (0)  in  a  given 

initial  state  i.   An  exponential  variate  with  mean  ■? — ; 

X±+M± 

is  generated  and  this  represents  the  length  of  time  that 
X„(t)  spends  in  state  i  before  jumping  to  state  i-1  or  i+1. 
In  order  to  decide  whether  XM(t)  next  visits  i-1  or  i+1,  a 

uniform  random  variate  p,  0  _<  p  <_   1 ,  is  generated.   If  p  is 

*i 
less  than  -r — ; ,  then  X,,(t)  visits  state  i+1  next  and  if  p 

Xi+^i        Mv 

i 
exceeds  -r—_ — 7   ^M(t)  goes  into  state  i-1.   Once  X„(t)  switches 

to  a  new  state  j,  we  repeat  the  process  of  determining  the 

Dength  of  X..(t)'s  stay  in  state  j  and  the  next  state  to  be 
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visited.   If  X.,(t)  is  in  state  0,  it  will  next  visit,  state 

M 

1  with  probability  one.   Similarly  if  X..(t)  is  in  state  M, 
it  will  next  visit  state  M-l  with  probability  one.   The  out- 
put of  this  program  will  be  discussed  together  with  that  of 
Program  I  in  the  next  section.   The  flow  chart  and  listing 
of  Program  II  are  included  in  Appendix  D. 

C.   DISCUSSION  OF  OUTPUT  OF  PROGRAMS  I  AND  II 

Both  simulation  programs  have  five  input  parameters, 
namely  the  total  number  of  machines  (N),  the  duration  of  each 
simulation  run  (TSTOP),  the  arrival  rate  of  messages  at  each 
machine  (A),  the  departure  rate  of  transmitted  messages  (u) 
and  the  initial  state  i.   Since  we  are  mainly  interested  in 
steady  state  results,  the  primary  output  of  both  simulation 
programs  is  the  probability  distribution  of  the  process 
XM(t)  when  in  steady  state. 

One  way  to  obtain  these  steady  state  probabilities  is  to 
let  the  simulation  program  run  for  t'  minutes  until  we  are 
reasonably  sure  that  steady  state  has  been  achieved,  and 

then  record  that  value  of  X„(t').   Since  each  simulation  run 

M 

produces  only  one  realization  of  XM(t'),  this  method  of  es- 
timating the  P.'s  proves  to  be  very  inefficient  because  of 
the  large  number  of  simulation  runs  required  to  produce  a 
consistent  estimate. 

A  better  method  to  estimate  the  P.'s  is  to  observe  the 

J 

time  average  probabilities  of  X.,(t)  when  it  is  already  in 
steady  state.  We  choose  a  starting  state  i  which  is  very 
close  to  the  steady  state  mean  My   so  that  XM(t)  will  reach 
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steady  state  shortly.   The  simulation  model  will  have  to  run 
for  a  considerable  amount  of  time  in  order  to  gather  enough 
sample  values.   However,  this  results  in  a  program  that  re- 
quires too  much  core  storage.   To  overcome  this  problem,  we 
made  200  simulation  runs  of  TSTOP  minutes  instead  of  a  single 
run  of  200  x  TSTOP  minutes. 

During  each  execution  of  the  simulation  program,  200 
simulation  runs  of  duration  TSTOP  minutes  each  are  made. 
Let  f .  denote  the  frequency  that  the  process  X„(t)  visits 
state  i.   At  the  beginning  of  the  program,  all  the  f.'s  are 
set  equal  to  zero.   During  each  simulation  run,  the  number 
of  messages  undergoing  transmission  is  observed  at  discrete 
one  minute  intervals.   If  this  number  is  k,  then  f.  is  incre- 
mented  by  one  and  all  other  f.'s,  j^k,  remain  unchanged. 

•J 

Thus,  during  the  entire  simulation  program,  X„(t)  is  observed 
at  200  *  TSTOP  points,  and  the  estimated  steady  state  proba- 
bilities are  given  by: 

f. 


^k    200xTSTOP 

Initially,  both  programs  are  tested  with  the  same  input 
parameters  to  ensure  that  they  are  yielding  results  that  are 
consistent  with  each  other.   From  then  on,  only  Program  II 
is  used  to  generate  the  various  output  that  are  used  to  val- 
idate the  theoretical  models,  because  time-wise,  it  is  much 
more  efficient  than  Program  I.   For  example,  with  N=50, 
A=2.5  and  u=10.0,  the  total  execution  time  of  Program  I  is 
about  18  times  that  of  Program  II. 
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VI.   SIMULATION  "NO  LOSS"  MODEL 

In  the  previous  sections,  we  have  developed  analytical 
and  simulation  models  for  a  system  of  N  machines  in  which 
messages  that  are  not  transmitted  immediately  are  lost.   A 
more  realistic  assumption  is  to  allow  these  messages  that 
are  not  transmitted  immediately  to  join  a  queue  at  the  re- 
spective machines.   We  call  a  model  based  on  this  less  strin- 
gent assumption  a  "no  loss"  model.   This  "no  loss"  model  is 
not  easily  amenable  to  theoretical  analysis,  but  it  can  be 
easily  simulated  by  a  program  written  in  GPSS .   The  assump- 
tions of  exponential  distribution  for  the  inter-arrival  and 
inter-departure  times  of  messages  could  now  be  relaxed  be- 
cause a  simulation  program  can  easily  handle  any  other  prob- 
ability distribution.   However,  the  previous  assumption  of 
exponential  inter-arrival  and  transmission  times  is  retained 
for  comparison  purposes. 

In  the  GPSS  program,  the  Sendox  machines  are  treated  as 
facilities  and  the  messages  that  arrive  at  the  various  ma- 
chines are  treated  as  transactions.   The  program  consists 
mainly  of  N  closely  similar  segments.   Each  segment  simulates 
the  activity  of  a  facility.   The  program  instructions  for 
each  segment  are  almost  identical  except  for  the  fact  that 
each  facility  has  a  different  facility  number  i,  i=l,2,...N. 

Arrivals  of  messages  are  simulated  through  the  generation 
of  random  arrivals  of  transactions  at  an  exponential  rate  A 
at  each  facility.   A  transaction  created  at  facility  i  is 
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randomly  assigned  to  a  destination  facility  j  ,  j^i,  inde- 
pendently of  other  transactions.   Then  facilities  i  and  j 
are  tested  to  see  whether  they  are  available  simultaneously. 
If  so,  the  transaction  proceeds  to  seize  both  facilities  for 
a  length  of  time  that  is  exponentially  distributed  with  rate 
u.   If  not,  the  transaction  will  have  to  wait  until  both 
facilities  i  and  j  are  simultaneously  free.   If  another  trans- 
action is  created  while  a  previous  transaction  is  either 
waiting  for  or  undergoing  service,  the  new  transaction  must 
queue  up.   Therefore  we  assume  a  first  come  first  serve 
queue  discipline.   This  procedure  is  followed  at  all  N  facil- 
ities. 

GPSS  is  preferred  to  FORTRAN  in  this  model  because  the 
logic  and  structure  of  a  GPSS  program  are  more  simple.   More- 
over queue  statistics  and  steady  state  probabilities  of 
XM(t)  can  be  gathered  with  little  extra  effort.   The  only 
shortcoming  is  that  this  program  cannot  adapt  easily  to 
changes  in  the  number  of  machines.   Thus  if  we  want  to  double 
the  number  of  machines  the  card  deck  will  almost  double  its 
size.   This  GPSS  program  was  only  written  for  N  equal  to  20, 
whereas  N  can  go  as  high  as  200  in  Programs  I  and  II.   A 
listing  of  this  program  is  included  in  Appendix  D. 
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VII.   COMPARISON  OF  RESULTS 

According  to  the  diffusion  "loss"  model,  when  X..(t)  is 
in  steady  state  X.,(t)  is  approximately  normally  distributed 
with  mean  x   and  variance  a2  given  by  equations  (4.16)  and 
(4.17).   Moreover,  using  x   and  a2  and  the  continuity  correc- 

o  S 

tion,  given  by  equations  (4.18),  (4.19)  and  (4.20),  the 
steady  state  probabilities  can  be  obtained  from  the  standard 
normal  table. 

We  are  mainly  interested  in  determining  the  accuracy  of 
the  diffusion  "loss"  model  for  various  values  of  N  and  p  in 
comparison  with  the  simulation  models.   For  ease  of  refer- 
ence equations  (4.16)  and  (4.17)  are  repeated  here. 

Xs  =  ^e  =  fye  f7'1' 

y|-ye<2-p)+1 

<   =   M   2(2+p-2y  )  (7-2) 


where   y   =  1  +  fp  -  /(l+ip)2-l 

Equations  (7.1)  and  (7.2)  implies  that  x   and  a2  depends 

only  on  the  ratio  of  X    and  u,  when  M  is  kept  constant.   To 

check  this  implication,  Simulation  Program  II  was  run  with 

five  sets  of  A  and  u  where  the  ratio  —  is  kept  equal  to  2. 

The  steady  state  means  and  variances  obtained  from  these 

runs  turned  out  to  be  extremely  close  to  x   and  o2  as  calcu- 

s       s 

lated  from  equations  (7.1)  and  (7.2)  and  these  results  are 
given  in  Table  I. 
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Table  I 


Steady  State  Mean  and  Variance  of  X„(t)  with 


p  constant,  N  =  100,  p  = 


_  y 


2A 


2.0. 


Arrival  rate 
of  messages 

X 

(#  per  hour) 

Departure  rate 
of  messages 

y 
(#  per  hour) 

Simulation 
"loss"  model 

Diffusion 
"loss"  model 

X 

s 

a2 
s 

X 

s 

a2 
s 

1.25 
2.5 

5.0 
10.0 
20.0 

5.0 
10.0 
20.0 
40.0 
80.0 

13.389 
13.406 
13.418 
13.410 
13.408 

7.712 
7.753 
7.787 
7.760 
7.760 

13.389 
13.  389 
13.389 
13.389 
13.389 

7.723 
7.723 
7.723 
7.723 
7.723 

(Each  simulation  runs  for  2400  hours) 
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Next,  keeping  M-10,  or  N=20,  six  more  runs  were  made  with 

the  values  of  p  ranging  from  0.5  to  10.   Again  x   and  a2 

s      s 

calculated  from  equations  (7.1)  and  (7.2)  were  in  close  agree- 
ment with  the  simulation  results.   These  results  are  summa- 
rized in  Table  II. 

Equations  (7.1)  and  (7.2)  also  imply  that  x   and  a2  are 

s       s 

linear  functions  of  N  keeping  p  constant  at  2.   Six  simula- 
tion runs  were  made  for  values  of  N  ranging  from  10  to  200. 
The  steady  state  probabilities  for  the  birth  and  death  "loss" 
model  were  also  calculated  for  the  same  values  of  p  and  N. 
Again,  the  results  indicated  that  all  three  models  agree  very 
well  as  far  as  steady  state  mean,  variance  and  cumulative 
probabilities  are  concerned.   The  steady  state  means  and 
variances  for  the  three  models  are  given  in  Tables  III  and 
IV  respectively.   In  Table  V  we  compared  the  cumulative 
probabilities  of  the  three  models  for  the  case  where  N=20 
and  p=2.   These  same  values  when  plotted  on  normal  probabil- 
ity graph  paper  (Figure  3),  gave  clear  visual  evidence  of 
their  close  agreement. 


33 


Table  II 


Steady  State  Mean  and  Variance  of  XM(t)  with 
N  constant,  N  =  40,  A  =  2.5. 


Departure  rate 

Simulation 

Diffusion 

of  messages 
y 

(#  per  hour) 

»-k 

"loss" 

model 

"loss" 

model 

X 

s 

a2 
s 

X 

s 

a2 
s 

2.5 

0.5 

9.976 

3.647 

10.000 

3.333 

5.0 

1.0 

7.631 

3.544 

7.600 

3.413 

7.5 

1.5 

6.263 

3.357 

6.277 

3.299 

10.0 

2.0 

5.352 

3.117 

5.359 

3.094 

30.0 

6.0 

2.537 

1.969 

2.540 

1.967 

50.0 

10.0 

1.676 

1.421 

1.678 

1.415 

(Each  simulation  runs  for  2400  hours) 
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Table  III 


Steady  State  Mean  of  XM(t)  with  X   and  y  constant, 
A  =  2.5,  y  =  10.0  and  p  =  ~   =  2.0. 


Number  of 
machines 

N 

St 

eady  State  Mean 

'   s 

Simulation 
"loss"  model 

Diffusion 
"loss"  model 

Birth  and  Death 
"loss"  model 

20 

2.703 

2.680 

2.697 

40 

5.373 

5.359 

5.376 

60 

8.049 

8.038 

8.055 

80 

10.722 

10.718 

10.735 

100 

13.397 

13.398 

13.414 

200 

26.803 

26.795 

26.811 

(Each  simulation  runs  for  2400  hours) 
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Table  IV 


Steady  State  Variance  of  X„(t)  with  X   and  u 
constant,  A  =  2.5,  u  =  10.0  and  p  =  —^-   =  2.0, 


Number  of 
machines 

N 

Steady  State  Variance  a2 

s 

Simulation 
"loss"  model 

Diffusion 
"loss"  model 

Birth  and  Death 
"loss"  model 

20 

40 
60 
80 

100 
200 

1.547 
3.052 
4.617 
6.142 
7.723 
15.440 

1.547 
3.094 
4.641 
6.188 
7.735 
15.470 

1.559 
3.106 
4.653 
6.199 
7.746 
15.484 

(Each  simulation  runs  for  2400  hours) 


36 


Table  V 

Steady  .State  Cumulative  Probabilities,  N=20, 
M=N/2=10,  X  =  2.5,  y  =  10.0,  p  =  ^  =2.0. 


State  j 

Steady  State  Cumulative 

Probabilities  P. 

.1 

Simulation 

Diffusion 

Birth  and  Death 

"loss"  model 

"loss"  model 

"loss"  model 

0 

0.0275 

0.0401 

0.0280 

1 

0.1646 

0.1716 

0.1681 

2 

0.4459 

0.4427 

0.4501 

3 

0.7485 

0.7453 

0.7470 

4 

0.9243 

0.9284 

0.9247 

5 

0.9865 

0.9887 

0.9864 

6 

0.9985 

0.9989 

0.9986 

7 

0.9999 

1.0000 

0.9999 

8 

1.0000 

1.0000 

1.0000 

9 

1.0000 

1.0000 

1.0000 

10 

1.0000 

1.0000 

1.0000 

(Each  simulation  runs  for  2400  hours) 
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Figure   3 
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VIII.   APPROXIMATION  OF  "NO  LOSS"  MODEL  BY  "LOSS"  MODEL 

Since  the  diffusion  "loss"  model  gives  relatively  simple 

formulas  for  x  ,  a2  and  P.'s,  it  appeared  desirable  to  re- 
s   s      j 

late  the  results  of  the  "loss"  model  to  the  "no  loss"  model. 
It  may  be  noted  that  in  the  "loss"  model,  only  a  fraction  of 
the  arriving  messages  are  transmitted,  whereas  in  the  GPSS 
"no  loss"  model  all  arriving  messages  are  eventually  trans- 
mitted. 

First,  we  need  to  determine  the  fraction  of  messages 
transmitted  in  order  to  obtain  the  effective  rate  of  message 
arrivals  at  a  "loss"  model.   In  the  "loss"  model,  the  prob- 
ability that  a  message  arrives  at  a  busy  machine  is  approxi- 
mately x   divided  by  M  or  y  ,  assuming  steady  state.   Since 
a  message  can  be  transmitted  only  if  both  the  transmitting 
and  receiving  machines  are  available,  the  probability  of  a 
message  being  transmitted  immediately  is  approximately 
(1-y  )2  assuming  independence  between  availability  of  ma- 
chines and  that  N  is  large.   Hence  X  „-  is  the  effective 

ef  f 

arrival  rate  of  messages  at  a  machine  in  the  "loss"  model  is 
given  by 

Xeff  =  Xl-V  C8.1) 

This  leads  to  the  idea  that  a  "no  loss"  model  with  input 
parameters  N,  X  and  u  can  be  approximated  by  a  "loss"  model 
with  input  parameters  N,  X'  and  u  where  X'  is  such  that 
X'(l-y  '  )2  =  X.   Using  equation  (4.15)  for  y  '  with  p'  =  £r-7 

G  G  ^j  A 
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the  equation  X'(l-y  '  )   =  X  can  be  rewritten  as 


X,[A*4XT)2-1  -  4X^2  =  X  <8'2> 

After  more  rearrangement  of  the  various  terms  (details  are 

given  in  Appendix  E),  we  obtain  an  equation  expressing  X' 

explicitly  as  a  function  of  u  and  X: 

X 
X'  =  —  ■  (8.3) 

(1-  i)2 

The  above  relationship  between  "loss"  and  "no  loss"  mod- 
els was  tested  using  the  simulation  models  with  five  differ- 
ent values  of  u  and  X1  as  input  parameters,  N  and  X  being 
kept  constant  at  20  and  2.5  respectively.   The  steady  state 
means  and  variances  are  summarized  in  Table  VI.   For  the 
case  of  N=20,  u=10.0,  X=2.5  and  X'=10.0,  the  steady  state 
cumulative  probabilities  are  also  calculated  using  the  dif- 
fusion "loss"  model,  and  the  results  from  all  three  models 
are  tabulated  in  Table  VII  and  plotted  in  Figure  4.   All  the 
results  indicate  that  the  approximation  of  the  "no  loss" 
model  by  the  "loss"  model  is  quite  accurate. 
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Table  VI 

Comparison  of  "no  loss"  model  with  its 
"equivalent"  "loss"  model,  N=20,  A=2.5, 


Simulation  "no  loss" 

model 

Simulation 

"loss  model 

X 

v 

X 

s 

a2 
s 

X' 

Xeff 

X 

s 

a2 
s 

2.50 

7.50 

6.677 

1.362 

22.500 

2.50 

6.694 

1.343 

2.50 

8.57 

5.758 

1.503 

14.407 

2.50 

5.862 

1.556 

2.50 

10.00 

4.972 

1.640 

10.000 

2.50 

5.031 

1.659 

2.50 

12.00 

4.105 

1.732 

7.347 

2.50 

4.192 

1.726 

2.50 

15.00 

3.319 

1.651 

5.625 

2.50 

3.357 

1.680 

(Each  simulation  runs  for  2400  hours) 
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Table  VII 

Steady  State  Cumulative  Probabilities  Calculated 
from  Three  Different  Models,  N=20,  y=10.0. 


Steady  State 

Cumulative  Probabilities  P. 

J 

State 
J 

Simulation 
"no  loss"  model 

Simulation 
"loss"  model 

Diffusion 
"loss"  model 

A=2.5 

X»=10.0 

X'=10.0 

0 

0.0000 

0.0001 

0.0002 

1 

0.0010 

0.0028 

0.0034 

2 

0.0240 

0.0245 

0.0266 

3 

0.1210 

0.1159 

0.1146 

4 

0.3480 

0.3347 

0.3492 

5 

0.659 

0.6380 

0.6507 

6 

0.8900 

0.8780 

0.8874 

7 

0.9810 

0.9790 

0.9734 

8 

0.9980 

0.9985 

0.9966 

9 

1.0000 

0.9999 

0.9998 

10 

1.0000 

1.0000 

1.0000 

(Each  simulation  runs  for  2400  hours) 
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IX.   CONCLUSION 

It  was  found  that  the  steady  state  results  of  the  diffu- 
sion and  the  birth  and  death  models  were  in  very  close  agree- 
ment with  the  simulation  model,  through  comparison  of  the 
steady  state  means,  variances  and  cumulative  probabilities. 
The  usefulness  of  the  diffusion  model,  as  compared  to  the 
birth  and  death  model,  lies  in  the  ease  with  which  its  steady 
state  mean  and  variance  and  its  discretized  probabilities 
can  be  computed  without  resorting  to  a  digital  computer. 

The  simple  form  of  equations  (7.1)  and  (7.2)  exhibits 
the  fact  that  the  steady  state  mean  and  variance  vary  line- 
arly with  N  and  depend  only  on  p,  and  not  on  the  individual 
values  of  X  and  u.   These  results  cannot  be  inferred  from 
the  steady  state  formulas  of  the  birth  and  death  model  or 
the  simulation  model. 

Lastly,  we  developed  a  simulation  model  for  the  much 
more  complicated  case  when  the  messages  are  allowed  to  wait 
for  transmission.   We  showed  that  a  "no  loss"  model  with  in- 
put parameter  N,  A  and  u  can  be  quite  accurately  approxi- 
mated by  a  "loss"  model  with  input  parameter  N,  X',  and  u 
where  X'  is  calculated  through  equation  (8.3)  from  the  other 
parameters. 
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APPENDIX  A 
Derivation  of  the  Characteristic  Function  of  S(t) 


Define  <f>„(9,t)  and  <J\,(8,t)  as  the  characteristic  func- 

M  M 

tions  of  X,,(t)  and  S  (t)  respectively: 

iex  (t) 

*M(6,t)  =  E[e    M    ]  .  (A.l) 

ies  (t) 

^M(6,t)  =  E[e        ]  (A. 2) 


Since 


XM(t)  =  My(t)  +  M   SM(t) 


we  obtain 

♦M(e,t)  =  e_9,/My(t)  ♦jjCe/ZH,  t)  (A. 3) 

Given  that  X  (t)  =  j,  equations  (4.3),  (4.4)  and  (4.5)  give 
the  respective  probabilities  that  X  (t+dt)  be  equal  to  j+1, 
j-1  or  j.  Hence  we  can  write  down  the  following  expression 
for  <!>M(e,t+dt): 

♦M(6,t+dt)  =  E[eieXM(t+dt)] 

=  E[eie(XM(t)  +  1){4M2-2M+2XM(t)(l-4M)+4XM(t)}|^I 

+  eieXM(t){l-2X,,(t)pXdt-(4M2-2M+2Xx,(t)(l-4M) 

+4X2(t))|^I}+  ei8(XM(t)-1){2XM(t)pAdt}  +  o(dt>] 

(A. 4) 
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Then  after  rearranging  the  terms  and  taking  the  limit  as 
t  -*■  °°,  we  obtain 

34>M(e.t)  =  lim  »M(e,t+dt)-«|»M(e,t)   • 

3t      dt+0        dt 

=  E[ei6XM(t)ei9{4M2-2M+2XM(t)(l-4M)+4X^(t)}2^ri 
+  eieXM(t){-2pAXM(t)-(4M2-2M+2XM(t)(l-4M)+4X^(t)) 
_}  +  ei0XM(t)e~i9{2pAXM(t)}]  (A. 5) 


2M- 
Further  rearranging  of  terms  results  in 


9^(6,t)  =  Erx2(t)eieXM(t)(eie-l)-^- 


+  2xxM(t)eiexMCt){iziM(eie_1)+p(e-ie_1)} 

+  2AMeieXM(t)(ei9-l)]  (A. 6) 


However,  since 

3<t>M(e,t) 


=  E[iX1I(t)ei8XM(t)] 


8  9        L   M 


32*M(e,t) 

§1F        "!  "M 


M,      =  E[-X2(t)ei9XM(t)] 


the  above  can  now  be  written  as 


a#M(e.t)   32*M(e,t)      4A 


ie 


(l-e^)^V  +  2XM<|>M(e,t)(e-  -1) 


3t         86'     v      y2M-l       Hi 


3<j>M(e,t)  4m-i  ie-  ,  ,  -ie   .        „. 

+  2\i ^ {2MTi(°  -l)-p(e   -1)1        (A. 7) 
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Rewriting   equation    (A. 3)    we   obtain   the   following   expressions 
for   $M(6 /fM, t )    and    its    first    and   second   partial    derivatives 

^(e/ZM.t)  =  ei0/i^(t\(e,t)  (A. 8) 

¥^j          ie'v%(t)   °V^j   +    .fli/B   ,,,,„    rfl   tU      ,.', 
^ =    e  *v     y{ ^ +    i8/My'(t)4'M(6,t)}       (A. 9) 


"M9^^  i9/My(t)    ^M(B,t) 

^o =    e  yv     y{ 5-5 +    i/My(t)^..(0,t)}  (A. 10) 


3< 

1  NT 


Mv    '        '     y  x6/My(t)r         Mv     '     '     ,     _  .  /tj    ,.  N       Mv     '     y 

8T2 =  e  { W1 +  2i/My(t)~ 53 

-My2(t)i|»M(e,t)}  (A. 11) 

If  we  rewrite  equation  (A. 7)  with  6//M  in  place  of  8,  we 
obtain 

9(frM(6/v/M,t)    _    8^M(9//M,t)  ie//jL4AM_ 

at  ae"2  (     e  J2M-i 

pl  A7.30M(9//"'t)/4M-l,     16//M    „        .    -16//M    _, 
+   2X/Mi g^ {2M=T(e  -1)-P(e         '       -1)} 

+    2AM4>M(e/»/M,t)(e:Le/v/^-l)  (A. 12) 

Using  equations  (A. 8)  through  (A. 11)  to  replace  all  the  terms 
involving  <f>  by  i> ,    and  factoring  out  e         ,  equation 
(A. 12)  becomes 
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N»M(e,t) 


^ +    i0/My'(t)^M(9>t) 

=    2M=l(1-e  ){ a?7 2i/My(t) My2(t)^M(e  ,t)} 

a  8 

oWM-riM-l,     ie//M1u    ,..       -ie//MNW8^M(9't) 
+   2A/Mi{-2Jjrj(e      '       -l)+p(l-e         '       )}{ — — 

+i/My(t)^M(6,t)}  +  2XM(el6//M-l)^M(e,t)  (A. 13) 

Substituting  e   '    and  e    '    with  their  power  series 

2 


i  -  ei9//M  =  zii  +  |i  +  0(1) 


2M    "\,3/ 


•M  M 


3/2 


1  _  e-i9//M  =  M  +  |i  +  0(-V) 
•M    2M      M3/2 

and  after  more  rearranging  of  terms  in  (A. 13)  we  may  write 
%*  *   19/BtM{y(t)-2XC^jy»(t)-(^  +p)y(t)-l)} 

-  (|fpj  -p)y(t)+l)2A(|^)^M+o(-4)  (A. 14) 

where  the  arguments  of  iK,(8,t)  have  been  dropped  for  ease  of 
exposition.   The  deterministic  process  y(t)  must  satisfy 
the  following  first  order  differential  equation,  as  we  ] et 

M  -*■  CO 

y!(t)  -2A  {y?(t)  -  (2+p)y(t)  -1}  =  0  (A. 15) 


48 


in  order  that  equation  (A. 14)  simplify  to 

||  =  -eff{2A(2+p-2y(t)}-  f%  {2X(y2(t)-(2-p)y(t)-l)} 

(A. 16) 

This  is  the  transformed  version  of  the  forward  differential 

equation  for  a  non-stationary  O.U.  process. 

We  have  to  solve  the  first  order  differential  equation 

(A. 15)  in  order  to  obtain  y(t)  as  explicit  function  of  t. 

For  the  sake  of  convenience,  y'(t)  and  y(t)  are  written  as 

dv 

-=*:  and  y  respectively.   Then  equation  (A.  15)  is 

||  =  2My2-y(2+p)+l}  (A. 17) 

Now  let 

v  =  1  +  ip 
and 

w2  =  (1+ip)2  -  1 

Then  equation  (A. 17)  can  be  written  as 

g  =  2A{(y-v)2-w2} 

which  in  turn  may  be  written  as 

( — ~-  )dy  =  4Awdt 

vy-v-w    y-v+wy  * 

Integrating  both  sides  of  the  equation  we  find 

In  y~V~W  =  4Xwt  +  Cj 
y-v+w 


or 


y-v+w       -4Xwt 

J =  c2e 

y-v-w 
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Assume  that  the  system  starts  off  empty,  hence  y=0  when  t=0, 
then  C2  is  evaluated  as 


c2 

= 

w-v 
w+v 

Now 

let 

Then 


w-v   -4Awt    y-v+w 

k  =  e      ==  

w+v  y-v-w 


1+k 

y  =  V  -        W 


Rewriting  y  as  a  function  of  t  again 

.,  ,  w-v  -4Xwt 
1  +  — — e 

y(t)  =  v  -  w  { ^^ ^Y-qr)  (A. 18) 

J v  '  H    w-v  -4Awt 

1  -   — 7— e 
w+v 

Further, 

limy(t)  /A-.r>\ 

ye  =  f"»  V_W  (A. 19) 
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APPENDIX  B 
Derivation  of  Stochastic  Differential  Equation 

From  equations  (4.21),  (4.22)  and  (4.23) 


dXM(t)  =  p  (x.t)dt  +  Z(t)/aM(x,t)dt  (B.l) 

{2M-2X  (t)}{2M-2X  (t)-l} 
gM(x,t)  =  M  2M-1 A  -2XM(t)pX      (B.2) 

{2M-2X  (t)}{2M-2X  (t)-l} 
aM(x,t)  =  2M_1 ~ A  +2XM(t)pX      (B.3) 


where  Z(t)  is  a  purely  random  Gaussian  process  with  zero 
mean  and  unit  variance.   Using  equation  (4.1)  we  can  express 
X„(t)  in  terms  of  S,,(t) 

M  M 

Xjj(t)  =  /MSM(t)  +  My(t) 
therefore 

dXM(t)  =  /MdSM(t)  +  My*(t)dt  (B.4) 

Using  equations  (B.2),  (B.3)  and  (B.4)  to  substitute  for 

3(x,t),  cc(x,t)  and  dX.,(t)  in  equation  (B.l)  we  obtain 

M 

/MdSM(t)  +  My' (t)dt 

{2M-2/MSM(t)-2My(t)}{2M-2v/MSM(t)-2My(t)-l} 
=  I 2M^1 Adt 

-  2pX{/MSM(t)+My(t)}dt]  +  Z(t)[2pX{/MSM(t)+My(t)}dt 

{2M-2/MS    (t)-2My(t)}{2M-2/MS    ( t )-2My ( t )-l}  , 

+  " 2MTI " Xdt]* 
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Carrying  out  all  the  necessary  multiplication,  rearranging 
some  of  the  terms  and  (since  we  are  eventually  going  to  let 
M  go  to  infinity)  approximating  2M-1  by  2M,  we  obtain  the 
following  expression  that  is  exact  in  the  limit. 

/MdS„(t)  +  My' (t)dt 
M 

1 

-  Xdt[2M-4My(t)+2My2(t)+SM(t){4/M(y(t)-l)+^|} 


+2SM2(t)-2p{/MSM(t)+My(t)}]  +  Z( t )/Xdt [ 2M-4My(t )+2My2 ( t ) 

1  i 

+SM(t){4/M(y(t)-l)+/M}+2SM2(t)+2p{/MSM(t)+My(t)}]2 

Dividing  by  /M  throughout  the  whole  equation  and  after  more 
rearranging  of  terms,  we  have 

dSM(t)    +    /My' (t)dt 

=    2/MAdt{y2(t)-(2+p)y(t)+l}+2XSM(t)dt{2y(t)-2-p+  |^} 

2X 
+  71   SM2(t)    +   Z(t)[2Adt{y2(t)-(2-p)y(t)+l} 

2XS    (t)dt  -  p,  , 

+  — 7i — {2^>-2+p+  Im>  +  IrVct)]*  <B-5> 

Now  we  let  M  go  to  infinity  and  assume  that  y'(t)  =  2X{y2(t) 
-(2+p)y(t)+l).   Then  with  S(t)  =  lim  S  (t)  we  have 

M->°o 

dS(t)  =  2X{2y(t)-2-p}S(t)dt+Z(t)[2X{y2(t)-(2-p)y(t)+l}dt]* 

Again  we  have  shown  that  the  process  S(t)  is  an  Ornstein- 
Uhlenbeck  process  with  the  following  infinitesimal  mean  and 
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variance,  conditioned  on  S(t)=s. 

3(s,t)  =  -3(t)s  =  -2A{2y(t)-2-p}s 

a(s,t)  =  a(t)  =  2A{y2(t)-y(t)(2-p)+l} 
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APPENDIX  C 
Derivation  of  Steady  State  Mean  and  Variance  of  X  (t) 


From  equation  (4.6)  we  have 

||  =  -3(t)eff  -  a(t)(|^H  (CD 

where 

3(t)  =  2X{2+p-2y(t)} 

ct(t)  =  2A{y2(t)-y(t)(2-p)+l} 

Let  the  mean  and  variance  of  S(t)  be  denoted  by  u(t)  and 
x(t),  and  since  S(t)  is  normally  distributed,  we  can  write 
its  characteristic  function  ip(8,t)  in  terms  of  u(t)  and 

T(t). 

^e^^J^^H^^t))  (c.2) 

Differentiating  ^(6,t)  with  respect  to  6  and  t  then  dropping 
the  arguments  of  ij;(  0  ,  t ) ,  we  obtain 

||  =  {iey-(t)-ie2T'(t)}e{ie^t)-4e2T(t)}        (C.3) 


at 


86 

w 


|$  =  {iy(t)-0T(t)}e{i^(t)-49  T^t)}  (C.4) 


e  substitute  these  equations  for  il>.  -^r ,    -r-^  into  equation 

do    ot 

(C.l).   Cancelling  out  the  factor  e< lQ^ ( t )-^9  x(t)}  we  ob_ 
tain 

ieu*(t)-£62T' (t)  =  -e(t)6{iy(t)-0T(t)}-ie2a(t)     (C.5) 


54 


By  equating  coefficients  of  iG  and  G 2  on  both  sides,  we  ob- 
tain the  following  two  first  order  differential  equations: 

u'(t)  =  -3(t)y(t)  (C.6) 

T'(t)  =  -2B(t)x(t)  +  a(t)  (C.7) 

Equation  (C.6)  can  be  easily  solved  by  separation  of  vari- 
ables [Ref.  7],  yielding 

t 
-/  3(x)dx 
y(t)  =  y(0)e  °  (C.8) 

It  is  reasonable  to  assume  that  the  noise  process  S(t)  starts 
off  at  time  zero  with  zero  mean  and  variance,  hence  y(0)=0 
and  t(0)=0. 

Substituting  the  initial  condition  y(0)=0  into  equation 
(C.8),  we  obtain 

y(t)  =  0       for  all  t  >_  0  (C.9) 

Equation  (C.7)  can  be  solved  [Ref.  7]  by  the  use  of  an  inte- 
grating factor,  giving  us  the  following  expression  for  x(t), 

with  t(0)=0 

t  z 

-2  /  6(x)dx   t      2   f    6(x)dx 

x(t)  =  e    °        {/  a(z)e   °        dz}        (CIO) 

o 

Note  that  8(t)  and  a(t)  are  functions  of  y(t)  where  y(t)  is 
a  complex  function  of  t,  given  by  equation  (A. 18): 


w_v      -4Awt 
1   +  — —  e 

y(t)    =   v   -   w    { ^±v __-}  (C.ll) 

J  v     '  \         w-v      -4A\vt 

1   -  — —  e 
w+v 
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Therefore  the  two  integrals  on  the  right  hand  side  of  equa- 
tion (CIO)  are  quite  intractable,  unless  t  is  very  large, 
in  which  case  we  are  able  to  derive  an  approximation  for 
T(t). 

For  the  typical  values  of  A  and  p  that  we  are  using,  the 

-4Awt 
term  e      in  equation  (C.ll)  vanishes  very  rapidly  with  t 

approaching  infinity.   For  example,  with  X=2.5  messages/hour, 
p=2  and  t-1.0  hour,  e"4Xwt  =  e^2'5**^1*  =  3.0  x  10"8. 
So  for  time  t  greater  than  one  hour,  y(t)  can  be  very  accu- 
rately approximated  by  y  where 


ye  =  l+|p  -  /(l+|p)2-l 

and  hence  3(t)  and  a(t)  can  be  approximated  by  the  following 
expressions. 

3(t)  -  2X(2+p-2ye)  =  Be  (C.12) 

«(t)  -  2A(ye2-ye(2-p)+l)  =  a&  (C.13) 

Referring  to  Figure  5,  we  can  see  that  when  z  is  very  large, 
such  that  the  area  under  the  curve  3(x)  up  to  time  x=1.0 
hour,  is  insignificant  when  compared  to  the  area  under  the 
curve  up  to  time  x-z ,  then  3  z  is  a  good  approximation  for 

z 

/  B(x)dx. 

0 
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X 

CQ 


> 

X 
«H 

o 

-p 
o 

ft 


m 

ho 

•H 
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z 

2  /  3(x)dx  23  z 

Similarly  a(z)e  can  be  approximated  by  a  e 

Thus  an  approximate  solution  for  r(t)  is 

a         -23  t 
x(t)  =  ■—-    [1  -  e   e  ]       t  >  10  hours 
23e 

-23  t 
However,  for  such  large  value  of  t,  e      ^  0,  so  we  have 

a 
t(t)  -  tS-  t  >  10  hours 

23e 

=  ye2-ye(2"p)+1 

2(2+p-2ye) 
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APPENDIX  D 
Notations  used  in  Simulation  Program  I 

a      arrival  time  of  next  message  at  machine  r 

S      departure  time  of  next  message  at  machine  r, 

S  =0  when  machine  r  is  not  in  use 
r 

L  denotes  the  machine  to  which  machine  r  is  transmitting 

N.  number  of  messages  undergoing  transmission  at  epoch  T. 

CLOCK  current  clock  time  of  simulation  model 

e  exponential  variate  with  mean  =  1 

X  arrival  rate  of  messages  at  each  machine 

u  departure  rate  of  a  message  undergoing  transmission 

N  total  number  of  machines 

n  discrete  uniform  (0,N)  variate 

TSTOP  length  of  one  simulation  run 
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r=l,2, . . .N 


)k 

ak=min{ar) 

S  .=min{S    ,S   ^0} 
j  r'    r7 


3L 


CLOCK=S . 
J 
%  ■  L. 
J 
S.=   0 
3 

V    ° 

N.=    N.-l 
l         l 

T.=    CLOCK 

l 
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Generate 
n,  n^k 


n=l,2, . . .N 


G> 


S.  =  CLOCK  +  e/y 
S  =  CLOCK  +  e/y 


N.=  N.+l 
1    1 

T.=  CLOCK 
l 


<S 


-© 


Output 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 


SIMULATION 
PROGRAM  I 


'LOSS*  MODEL 


SIMULATE  A  SYSTEM  OF  N  SENDOX  MACHINES 
ASSUMING  THAT  ANY  ARRIVAL  THAT  FAILS 
TO  GET  ON  THE  SYSTEM  IS  LOST 


A(J) 
S(JJ 

LINKt J) 


NSYS 

N 
M 

TIME 

DELT 

P4 

P5 

PCF4 
PDF5 
ARPM 

DEPM 

TSTOP 
NRUN 

IX 

CLOCK 

NARR 

NCN 


TIME 

TIME 

=  0 

CONT 

IS  R 

TO  M 

NUMB 

TOTA 

MAX  I 

BE  T 

RECO 

DEPA 

DISC 

RELA 

TIME 

DISC 

CUMU 

CUMU 

RECI 

EACH 

RECI 

TRAN 

LENG 

NUMB 

RAND 

CURR 

NUMB 

NUMB 

TRAN 


OF 

OF 
IF 
AIN 
ECE 
ACH 
ER 
L  N 
MUM 
RAN 
RDS 
RTU 
PET 
TIV 

AV 
RET 
LAT 
LAT 


PRO 

SMI 

TH 

ER 

CM 

ENT 

ER 

ER 

SMI 


NEXT 

NEXT 
MACHIN 

THE  N 
I  VING 
I  HE  J 
OF  BUS 
UM6ER 

NUM3E 
SMITTE 

THE  I 
RE 

E  TIME 
E  FREQ 
ERAGE 
E  TIME 
IVE  PR 
IVE  PR 
CAL  OF 
CHINE 
CAL  OF 
TTED  M 
OF  EAC 
OF  SIM 
NUMBER 

CLOCK 
OF  MES 
OF  MES 
TTED 


ARRIVAL  AT  MACHINE  J 
DEPARTURE  AT  MACHINE  J 
E  IS  NOT  IN  USE 
UMBER  OF  THE -MACHINE  THAT 
A  MESSAGE  OR  TRANSMITTING 

Y  MACHINES 

OF  MACHINES 

R  OF  MESSAGES  THAT  CAN 

D  ,  =  N/2 

NSTANT  OF  ARRIVAL 


INT 
UENC 
PROB 

PRO 
OBAB 
08AB 

ARR 
IN  U 

DEP 
ESSA 
h  SI 
ULAT 

SEE 

TIM 
SAGE 
SAGE 


ER 

Y 

A6 

BA 

IL 

IL 

IV 

NI 

AR 

GE 

MU 

ID 

D 

E 

S 

S 


OR 


VAL  AT  WHICH 

ARE  COLLECTED 

ILITY 

BILITY 

ITY  FOR 

ITY  FOR 

AL  PATE 

TS  OF  MINUTES 

TURE  RATE  CF 

IN  UNITS  OF  MINUTES 
LATIuN  RUN  IN  MINUTES 
N  RUNS 


P4 
P5 

OF 


MESSAGE  AT 


THAT 
THAT 


ARRIVE 
ARRIVE 


AND    ARE 


C 
C 
C 


DIMENSION    A(200),     S(200),     LINK(200) 
2t    KSYS(5000),     TIME(5000) 

3,  PK101),     P2QG1),     P3Q01),    Y(101) 

4,  DELT(IOOO) 

5,  P4Q01),     P5(101) 

6,  PDF4C101J,     PDF5Q01) 
CALL    OVFLOW 

INITIALIZATION 

READ  (5,1000)  NRUN 

WRITE(6,1000J  NRUN 

READ   (5,1000)  N 

WRITE     (6,1000) 
1000    FC P.MAT     (15) 

READ       (5,2000) 

WRITE     (6,2000) 
2000    FORMAT     (3F10.2) 

Ml=N/2+l 

DC     10    1=1, Ml 

P4( I )=0.0 

P5( I )=0.0 
10      CONTINUE 

IX=52841 

START    OF    SIMULATION    »UNS 


N 

ARRM, DEPM, TSTOP 
ARRM, DEPM, TSTOP 
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c 
c 
c 
c 
c 


c 
c 
c 
c 


c 
c 
c 
c 


20 


25 
30 


C 
C 
C 


40 


C 
C 
C 
C 


50 


DO  95  IC=ltNRUN 
IX=IX+IC*246 

CLOCK=0.0 

NARR=0 

NCN=0 

JAD  =  1 

NSYS(JAD)=0 

TIME( JAD)=CLOCK 


AT  THE  START 
AND  SET  A(K) 


1,N 


OF  EACH  SIMULATION  RUN   SET  ALL  S(K)=0 
EQUAL  TO  TIME  OF  FIRT  ARRIVAL 


DC  20  I 

S(I)=0 

CALL  EXPON  (IX,E,1) 

Ad  )=E*ARRM 

CONTINUE 

DO  25  1=1, Ml 

P1CI)=0.0 

P2(I)=0.0 

CONTINUE 

CONTINUE 


SUBROUTINE 
SUBROUTINE 


ARNEXT 
ARNEXT 


SEARCHES 
SEARCHES 


CALL  ARNEXT  (A,K,N) 
CALL  DPNEXT  (S,L,N,SMJ 
IF  (A(K)  .GE.  SM)  GO  TO  50 
IF  (A(K)  .GE.  TSTOP)  GO  TO 


FOR 
FOR 


60 


MIN  A(K) , K  =  l  ,N 

MIN  S(K) ,K=1,N,S(K) tfO 


THE  NEXT  EVENT  IS  AN  ARRIVAL  AT  MACHINE  K 
TEST  TO  SEE  WHETHER  IT  WILL  BE  TRANSMITTED 

CLCCK=A(K) 

CALL  EXPON  < IX, E, 1  ) 

A(K)=A(K)+E*ARRM 

NARR=NARR+1 

IF  (S(K)  .GE.  CLOCK) 

CONTINUE 


GO  TO  30 


GENERATE  DESTINATION  MACHINE  NDES=1 , N , NDES #K 

CALL  RANDOM  (IX,Z,1) 

NDES=Z-N 

NDES=NDES+1 

IF  (NDES  .EQ.  K)  GO  TO  40 

IF  (S(NDES)  .GE.  CLOCK)  GO 


CALL  EXPON  (IX,E,1) 

E=E*DEPM 

S(NDES)=CLOCK+E 

LINK(NDES)=K 

S(K)=CLOCK+E 

LINK(K) =NDES 

JAD=JA0+1 

NCN=N0N+1 

NSYS( JAD)=NSYS( JAD-1J+1 

TIME( JAD)=CLOCK 

GO  TO  30 

CONTINUE 


TO  30 


THE  NEXT  EVENT 
AND  MACHINE  JL 


IS  A  DEPARTURE 
(=LINK(D) 


AT  MACHINE  L 


IF  (SM    .GE.  TSTOP)  GO  TO  60 

JL=LINK(L ) 

S(JL)=0 

CLOCK=S(L) 

S(L)=0 

JAL=JAD+1 

NSYS(  JAD)=NSYS(  JAD-D-1 
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60 


C 
C 
C 


C 
C 
C 


70 


75 


80 


90 
95 


C 
C 
C 
C 


96 
97 


8200 
99 

8800 


C 
C 
C 


TIME< JAD)=CLOCK 

GC    TO    30 

CONTINUE 

TIME( JAD+l)=TSTOP 

COMPUTE    TIME    AVERAGE    PROBABILITIES 

DO    70    1=1, JAD 

I  I  =  NSYS (I)+l 

P1(II)=P1(II)+TIME(H-1)-TIME(I) 

CONTINUE 

COMPUTE    DISCRETE    TIME    PROBABILITIES 

JJ  =  1 

NSTGP=TSTOP+0. 1 

DO  30  I=1,NSTCP 

DELT( I )=FLOAT(  I  ) 

CONTINUE 

JJ=JJ+1 

IF  (DELT(I)  .GT.  TIME(JJ))  GO  to  75 

JJ=JJ-1 

KK=NSYS(JJ )+l 

P2(KK)=P2(KK)+1.0 

CONTINUE 

DO  90  1=1, Ml 

Yd  )  =  FL0AT(I-1) 

P1(I)=P1(I J/TSTQP 

P2(I)=P2( I )/TSTOP 

P4(I)=P4( I )+Pl( I ) 

P5(I)=P5(I  )  +  P2( I  ) 

CONTINUE 

IX=52841 

CONTINUE 

SP4=0.0 

SP5=0.C 

VP4=0.0 

VP5=0.0 

COMPUTE  STEADY  STATE  MEAN,  VARIANCE  AND 
CUMULATIVE  PROBABILITIES  USING  P4  L    P5 

DO  99  1=1, Ml 
P4(I)=P4( I )/NRUN 
P5(IJ=P5( I j/NRUN 
SP4=SP4+P4(I )*( 1-1) 
SP5=SP5+P5(I 
VP4=VP4+P4(I 
VP5=VP5+P5(I 
IF  (I  .NE.  1) 
PDF4(I)=P4(I ) 
PDF5( I )=P5(I ) 
GO  TO  97 
CONTINUE 

PDF4U  )  =PDF4(I-1)+P4(  I  ) 
PDF5( I  )=PDF5( I-l)  +  P5( I ) 
CONTINUE 
I  1  =  1-1 

WRITE  (6,8200)  II 
1,  PDF5(I) 
FCRMAT  lI10,5F12-« 
CONTINUE 
WRITE  (6,3800) 
FORMAT  (////) 
VP4=VP^-SP4*#2 
VP5=VP5-SP5**2 


)  *  (  I  - 1  ) 
)*(I-1)*(I-1) 

)*( I-1)*(I-1) 


GO  TO  96 


'+, 


P3(  I)  ,  P4( I)  ,  P5( I ) ,  PDF4( I) 
/) 


7000 


PRINTOUT  OF  RESULTS 

WRITE(6,7000)  NRUN,SP4,VP4 

FORMAT  llOXt'TIME  AVERAGE  PROBABILITY 


1  '  SIMULATION  RUNS',  /  ,  10X,  •  MEAN 


WITH' 
'  ,F10.4,  ' 


,15, 
VARIANCE  =' 
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2  F10.4,//) 

WRITE(6,3000)  NRUN,SP5,VP5 
8000  FORMAT  (10X,  'DISCRETE  TIME  PROBABILITY  WITH', 15. 

1  '  SIMULATION  RUiSiS't/t  10X,  'MEAN  =  «,F10.4, 

2  •   VARIANCE  =» , F 10.4, //// ) 
CALL  PLOTP  (YtP3,Ml,0J 
WRITE  (6,8800) 

CALL  PLOTP  <Y,P4,M1,0) 

WRITE  (6,8800) 

CALL  PLOTP  (Y,P5,M1,01 

WRITE  (6,8800) 

STOP 

END 

SLBROUTINE  ARNEXT  (A,K,N) 

DIMENSION  A(100) 

K  =  l 

AM=100000.0 

DO  50  1=1, N 

IF  (AC  I  )  .GE.  AM)  GO  TO  50 

AM=A(I) 

K=I 
50   CONTINUE 

RETURN 

END 

SUBROUTINE  DPNEXT  (S,L,N,SM) 

DIMENSION  S(IOO) 

L  =  l 

SM=2000 000.0 

DO  50  1=1, N 

IF  (S(I)  .LE.  0.00001)  GO  TO  50 

IF  (S( I )  .GE.  SM)  GO  TO  50 

SM=S(I ) 

L  =  I 
50   CONTINUE 

RETURN 

END 
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Notations  used  in  Simulation  Program  II 

a  arrival  rate  of  messages  at  each  machine 

d  departure  rate  of  a  message  undergoing  transmission 

N  total  number  of  machines 

M  =  N/2 

X.  rate  of  entering  state  i 

u .  rate  of  leaving  state  i 

1ST  initial  state  at  time  0 

N.  number  of  messages  undergoing  transmission  at  epoch  T. 

e  exponential  variate  with  mean  =  1 

q  uniform  (0,1)  variate 

TSTOP  length  of  one  simulation  run 
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Ak 


T.      =T.+eA. 
x+1      1  1 


© 


Read 
N,  a,  d 

1ST 


AZ 


M   =   N/2 
Ml    =   M+l 


AIL 


Calculate 
X . ,    i=l , . . . M 

Vy    j=2,. ..Ml 


ALL 


1   =    1 

T.    =    0 

l 

N.    =    1ST 

l 


)k 

Ti+l=VeAM 
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© 


e 


JsL 


i   =    i+1 

N.=N.    n+l 
1      l-l 


\k 


Xi+1      i      X.+y. 


Q^^d 


© 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 


SIMULATION 
PROGRAM  II 


•LOSS'  MODEL 


TO  SIMULATE  THE  TRANSITION  OF  STATES  IN 
BIRTH  AND  DEATH  PROCESS 


NRUN 

M 

RARP 

RDEP 

RL(J) 
RM{J1 

IX 

TSTOP 

NSYS(J) 

1ST 
Pit  P4 

P2,  P5 

PDF4 
PDF5 


NUMB 

MAXI 

CAN 

RATE 

MACH 

RATE 

UNDE 

RATE 

RATE 

RAND 

LENG 

NUMB 

TRAN 

INIT 

STEA 

FOR 

STEA 

FOR 

CUMU 

CUMU 


ER  OF  S 
MUM  NUM 
BE  TRmN 

OF  APR 
INE 

OF  DEP 
RGOING 

OF  ENT 

OF  LEA 
OM  NUMB 
TH  of  e 
ER  OF  M 
SMISSIO 
IAL  STA 
DY  STAT 
li  200 
DY  STAT 
li  200 
LATIVE 
LATIVE 


IMULATICN  RUNS 

3ER  OF  MESSAGES  THAT 

SMITTED 

IVAL  OF  MESSAGES  AT  EACH 

ARTURE  OF  A  MESSAGE 

TRANSMISSION 

ERING  STATE  J 

VING  STATEJ 

ER  SEED 

ACH  SIMULATION  RUN  IN  MINUTES 

ESSAGES  UNDERGOING 

N  AT  EPOCH  T(J) 

RTING  STATE 

E  TIME  AVERAGE  PROBABILITIES 

SIMULATION  RUNS  RESPECTIVELY 

E  DISCRETE  TIME  PRC3 ABI L I T I ES 

SIMULATION  RUNS  RESPECTIVELY 

STEADY  STATE  PROBABILITIES  P4 

STEADY  STATE  PROBABILITIES  P5 


C 
C 
C 

c 


DIMENSION  RL(lOl),  RM(lOl) 
If  NSYS(5000),  T(5000) 

2,  PK101),    P2(101),    P3(101)f 

3,  DELT(IOOO) 

4,  PDF4C101),    PDF5(101) 

5,  P4(101)f    P5(101) 
CALL    OVFLOW 

INITIALIZATION 


READ  (5,1000)  NRUN 

WRITE(6,1000)  NRUN 

READ  (5,1000)  M 

WRITE(6,1000)  M 

READ  (5,1000)  1ST 

WRITE(6,1000)  1ST 
00  FORMAT  (15) 

READ  (5,2000)  RARR,  RDEP,  TSTOP 

WRITE<6  ,2000)  RARR,  RDEP,  TSTOP 
00  FGPMAT  (3F10.4) 

RARR=1.0/RARR 

RDfcP=l .O/RDEP 

M1=M+1 

M2=2*M 

RARR=RARP/(M2-1 ) 

P4(l)=0.0 

P5(l)=0.0 

COMPUTE  RATES  OF  ENTERING  AND 
LEAVING  THE  VARIOUS  STATES 


DC  10  1=2, Ml 

I2=2*( 1-2) 

RL( 1-1) =(M2-I2)*RARR*(M2-I2-1J 

RM( I)  =  (  I-1)#RDEP 
P4( I )=0.0 
P5(I)=0.0 
10   CONTINUE 
IX=35641 


X(101) 


10 
20 
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c 
c 

c 


c 
c 
c 
c 


c 
c 
c 
c 


c 
c 
c 
c 


c 
c 
c 


c 
c 
c 


c 
c 
c 


c 
c 
c 


15 


20 


30 


40 


50 


60 


70 


START    OF    SIMULATION    RUNS 
REPEATED    NRUN    TIMES 

DO    95    NR=lfNRUN 
IX=IX+NR*513 
DO    15    1=1, Ml 
Pl(  I)=0.0 
P2(I)=0.0 
CONTINUE 

AT  THE  BEGINNING  OF  EACH  RUN  SET  TIME  T(0)=0 
AND  NSYS(0)=IST 

IND  =  1 

T(IND)=0.0 
NSYS( IND)=1ST 

CONTINUE 

CALL  EXPON  { IX, E, 1) 

J=NSYSl IND)+1 

IF  (J  .EQ.  1)  GO  TO  30 

IF  (J  .EQ.  Ml)  GO  TO  40 

RATE=RL( J)+RM( J ) 

T(IND+1)  DENOTES  THE  INSTANT  WHEN  THE  PROCESS 
CHANGES  STATE 

TdND  +  1  )  =  T(IND)+E/RATE 

IF  (T(IND+1)  .GE.  TSTCP)  GO  TO  70 

TO  DETERMINE  WHETHER  PROCESS  GOES  TO  STATE 
J+l  OR  J-l 

PB=RL< J )/RATE 

CALL  RANDOM  (  I  X  ,  Y  ,  1 ) 

IF  (Y  . LE.  PBJ  oO  TO  50 

GO  TO  60 

CONTINUE 


PRGCESS  IS  IN  STATE  0  ALWAYS  GOES  TO  STATE  1 

70 


TdND  +  1  )=T(IND)+E/RL( J) 

IF  (TdND+1)  .GE.  TSTOP)  GO  TO 

GO  TO  50 

CONTINUE 


PROCESS  IS  IN  STATE  M  ALWAYS  GOES  TO  STATE  M-l 

TO  70 


TdND+1  )=T(  IND)+E/RM(  J  ) 

IF  (T(IND+1)  .GE.  TSTOP)  GO 

GO  TO  60 

CONTINUE 


PROCESS  GOES  TO  STATE  J+l 

K=NSYS( IND)+1 
IN0=IND+1 

NSYS( IND)=K 

P1(K)=P1(K)+T(IND)-T( IND-1) 
GC  TO  20 

CONTINUE 

PROCESS  GOES  TO  STATE  J-l 

K=NSYS( INDJ-1 
IND=IND+1 

NSYS( IND)=K 

PliK+2) =Pl(K+2)+T{ INU)-T( IND-1) 

GC  TO  20 

CONTINUE 

K  =  NSYS(  IND)  +  1 

PHK)=P1(K)  +  TST0P-T(IND) 
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C      COMPUTE  P2  THE  DISCRETE  TIME  PROBABILITIES 
C 

NSTOP  =  TSTOP+0.  1 

JK=1 

DC  80  I =1, NSTOP 

DELTd  )  =1 
75   CONTINUE 

JK=JK+1 

IF  (DELT(I)  .GE.  T(JK))  GO  TO  75 

JK=JK-1 

KK=NSYS(JK)+1 

P2(KK)=P2<KK)+1 
80   CONTINUE 

DO  90  1=1,  Ml 

X(I )=I-1 

P1(I)=P1( I )/TSTOP 

P2(I)=P2< I )/TSTOP 

P4(I)=P4( I  )+Pl(T  ) 

P5(I)=P5( I )+P2(  I  ) 
90   CONTINUE 

IX=35641 

95  CONTINUE 
SP4=0.0 
SP5=0.0 
VP4=0.0 
VP5=0.C 

C 

C      CALCULATE  STEADY  STATE  MEAN,  VARIANCE  AND 

C      CUMULATIVE  PROBABILITIES 

C 

DC  99  1=1, Ml 

P4(  I)  =  P4< I J/NRUN 

P5(I)=P5< I J/NRUN 

SP4  =  SP<+  +  P4(I  )*(  1-1) 

SP5=SP5+P5 (I )*( 1-1) 

VP4=VP4+P4(I l*LI— 1 J*C I-ll 

VP5=VP5+P5(I )*(I-1)*U-1) 

IF  (I  .NE.  1)  GO  TO  96 

PDF4U)  =  P4(I  ) 

PDF5( I )=P5(I) 

GO  TO  9  7 

96  CONTINUE 

PDF4( I) =P0F4( I-l)+P4( I ) 
PCF5( I )=PDF5( 1-1 J+P5( I ) 

97  CONTINUE 
11=1-1 

WRITE  (6,8200)  II,  P3  (  I  ) ,  P4U),  P5(I),  PDF41I) 
1,  PDF5(I) 
8200  FORMAT  ( I  10, 5F12 . 4 , / ) 
99   CONTINUE 
C 
C      PRINTOUT  OF  RESULTS 


C 


WRITE  (6,8800) 
8800  FORMAT  (////) 

VP4  =  VP^*-SP4**2 

VP5=VP5-SP5**2 

WRITE(6, 7000J  NRUN,SP4,VP4 
7000  FORMAT  (10X,'TIME  AVERAGE  PROBABILITY  WITH', 15, 

1  '  SIMULATION  RUNS' ,/,  10X, 'MEAN  =  '  , F 10 . 4 , ' VAR I ANCE  =  ' 

2  F10.4,////) 
NRITE(6,8000)  NRUN,SP5,VP5 

8000  FORMAT  (10X,  'DISCRETE  TIME  PROBABILITY  WITH', 15, 

1  •  SIMULATION  RUNS' ,/ ,10X, 'McAN  =  «,F10.4, 

2  '   VARIANCE  = ' , F 1 0.4 , //// ) 
CALL  PLOTP  (X,P3,M1,0) 
WRITE  (6,8600) 

CALL  PLCTP  (X,P4,M1,0) 

WRITE  (6,8800) 

CALL  PLOrp  IX,P5,M1,0) 

STOP 

END 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


THIS  PROGRAM  CALCULATES  THE  STEADY  STATE  MEAN , VAR I ANCE 
AND  CUMULATIVE  PROBABILITIES  OF  A  BIRTH  AND  DEATH 
PROCESS 


M 
P 

CL 

CM 

PDF 

XMEAN 

XVAR 


MAXIMUM  NUMBER  OF  STATES 

VECTOR  OF  STEADY  STATE  PROBABILITIES 

ARRIVAL  RATE  OF  MESSAGES  AT  EACH  MACHINE 

DEPARTURE  RATE  OF  A  MESSAGE  WHILE 

UNDERGOING  TRANSMISSION 

CUMULATIVE  STEADY  STATE  PROBABILITIES 

STEADY  STATE  MEAN 

STEADY  STATE  VARIANCE 


DIMENSI 
READ(5, 

1000  FCRMAT 
DC  10  K 
READ  (5 
WRITE (6 
P(l)=l. 
SUM  =P( 
M1  =  M+1 
M2=2*M 
CL=CL/( 
DC  20  I 
J2  =  2=M I 
C1=(M2- 
C2=(  1-1 
P(  I)  =  P( 
SLM=SUM 
20  CONTINU 
XMEAN=0 
XVAR=0. 
WRITE  ( 

2000  FCRMAT 
DO  60  I 
P(I)=P( 
X ( I )=I- 
XMEAN=X 
XVAR=XV 
IF  (I  . 
PDF(I)- 
GO  TO  4 
30  CCNTINU 
PCF(I  )  = 
40  CONTINU 
11=1-1 
WRITE  ( 

3000  FCRMAT 

60   CCNTINU 

XVAR=XV 

WRITE  { 

WRITE  ( 

4000  FCRMAT 
1,  F10.4 
CALL  PL 
WRITE  ( 
10  CONTINU 
STOP 
END 


ON  X(201) ,  P(201), 
1000)  NCARD 
(I5,2F10.4) 


PDF(201) 


=1, NCARD 
,1000)  M, 


,1000) 

0 

1) 


M, 


CL,  CM 
CL,CM 


M2-1) 

=  2,  Ml 

-2) 

J2 J*CL*(M2-J2-1) 

)*CM 

I-1)*C1/C2 

+  P(I) 

E 


0 

6, 
(/ 

=  1 
I) 
1 

ME 
AR 
NE 
P( 
0 


2000) 
////) 
i  Ml 

/SUM 

AN+Pd  )*(  1-1  ) 
+  P(I)=M  I-l)*(  I- 

.  1)  GO  TO  30 
I  ) 


1) 


PDF(I-1)+P( I ) 
E 


6,3 
(II 
E 

AP- 
6,2 


.  s 


(10 
,// 

OTP 
6,2 
E 


000)  II,  P(  I) ,  PDF(I ) 
0,2F15.4,/) 

XMEAN**2 

000) 

QUO)  XMEAN, XVAR 

Xt'MEAN  =  ■  ,F10. 4, 10X,  'VARIANCE 

//) 

(X,P,M1,0) 
000) 


72 


Listing  of  GPSS  Program 

This  program  simulates  the  arrivals  and  departures  of 
messages  at  a  system  of  20  sendox  machines.   Functions  EXP 
generates  exponential  variates  with  mean  of  one  and  UNI  gen- 
erates discrete  uniform  variates  between  1  and  20.   The  pro- 
gram consists  of  20  closely  similar  segments.   Each  segment 
simulates  the  activities  at  one  machine.   At  the  end  of  each 
simulation  run,  this  program  gives  the  expected  waiting 
times  at  each  machine  and  steady  state  probabilities  of 
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APPENDIX  E 
Calculation  of  A' 


In  the  previous  sections,  we  had  shown  that  a  "no  loss" 

model  with  input  parameters  N,  A  and  u  can  be  approximated 

by  a  "loss"  model  with  input  parameters  N,  A'  and  u,  where 

A1  is  chosen  such  that  the  effective  arrival  rate  A'(l-y')2 

e 

for  the  "loss"  model,  is  equal  to  A. 

When  X.,(t)  is  in  steady  state,  the  mean  number  of  busy 
machines  is  Ny  ,  hence  the  probability  of  finding  a  busy 
machine  is  Ny  /N  or  y  .  In  the  "loss"  model,  a  message  is 
transmitted  only  if  it  finds  both  the  transmitting  and  re- 
ceiving machines  free.  The  probability  of  this  event  is 
approximately  (1-y  )2  assuming  independence  between  avail- 
ability of  machines  and  a  large  N.   Therefore,  the  effective 

arrival  rate,  A  ._.  for  a  "loss"  model  is  A '(1-y')2  where 
'   ef  f  e 

y'  is  given  by  equation  (4.15)  as 


y'  =  1+ip'  -  /(l+ip')2-l 


with  p'  =   i: 


2A' 
Therefore,  we  require  that 


A'(/(l+ip'  )2-l  -  hQ'  )2  «  A  (E.l) 

After  expanding  the  quadratic  term  in  equation  (E.l)  and  re- 
arranging terms,  we  obtain 

A 


(1+iP')  -  t47T  =  /p'(l  +  iP'  )  (E-2) 

A   p 
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Squaring  both  sides  of  equation  (E.2)  and  recalling  that 
p=  — ,  we  obtain 

p'  =  p  +  -  -  2 

.       P 

Making  the  substitution  for  p  and  p',  we  finally  have 


A'  =  — 

(l-  V 


Using  this  rate  of  arrivals  in  our  "loss"  model  we  are  able 
to  approximate  the  system  size  characteristics  of  the  "no 
loss"  model  with  arrival  rate  X. 
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